rm(list = ls())

wd <- ".../Replication/"
setwd(wd)

# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# Load data ----
d <- read.csv("Data/baseline database.csv")
data <- subset(d, d$group!="Liberado por eleccion")


ITT_mechs <- list()

ITT_mechs[['Robbery']] <- with(data, lm_robust(robbery~lottery_winner))
ITT_mechs[['Assault']] <- with(data, lm_robust(assault~lottery_winner))
ITT_mechs[['Discretionary']] <- with(data, lm_robust(discretionary~lottery_winner))

modelsummary(ITT_mechs, stars = T, output = "latex")

declare <- declare_ra(N = nrow(data), m=sum(data$lottery_winner)) 
ri_mechs <- list()
ri_mechs[['Robbery']] <- summary(conduct_ri(formula = robbery ~ lottery_winner,
                                            declaration = declare,sharp_hypothesis = 0,
                                            assignment = "lottery_winner",
                                            data = data,sims = 10000
))


ri_mechs[['Assault']] <- summary(conduct_ri(
  formula = assault ~ lottery_winner,
  declaration = declare,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = data,
  sims = 10000
))


ri_mechs[['Discretionary']] <- summary(conduct_ri(
  formula = discretionary ~ lottery_winner,
  declaration = declare,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = data,
  sims = 10000
))


mylabels <- c(paste0("Enslaved\n(N=", sum(data$lottery_winner==0), ")"), 
              paste0("Emancipated\n(N=", sum(data$lottery_winner==1), ")"))
interval1 <- -qnorm((1-0.90)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# a. Robbery ----

figdata <- Rmisc::summarySE(data, 
                            measurevar="robbery", 
                            groupvars=c("lottery_winner"))

figdata$point <- sprintf(figdata$robbery*100, fmt = '%#.2f')

# ri p-val
ITT_CACE_pri <- paste0(#"CACE: ", round(CACE_mechs[['Robbery']]$coefficients[2],3)*100,"%p\n",
  "ITT: ", round(ri_mechs[['Robbery']]$estimate,3)*100,"%p\n",
  "(p=", round(ri_mechs[['Robbery']]$two_tailed_p_value,3),")")



robbery <- ggplot(figdata, aes(x=factor(lottery_winner), 
                               y=100*robbery, 
                               fill=factor(lottery_winner))) +
  geom_bar(stat="identity", width=0.45, alpha=0.9) +
  geom_text(aes(label = point), size = 4.0, vjust=-0.5, hjust = 0.4, nudge_x = 0) +
  xlab("") +
  ylab("Incarcerated (%)") +
  ylim(-2, 27) +
  scale_x_discrete(breaks=c("0", "1"), labels=mylabels)+
  theme_minimal() +
  theme(axis.text.y=element_text(colour = "gray30", size=11), axis.ticks.y=element_blank()) +
  theme(axis.text.x=element_text(colour = "gray30", size=11), axis.ticks.x=element_line(colour="gray30")) + 
  geom_segment(aes(x=1, y=19, xend=1.2, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.2, y=21, xend=1.8, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.8, y=21, xend=2, yend=19), colour="darkgray") + 
  annotate("text", x=1.5, y=21, label=ITT_CACE_pri, size=4) +  
  scale_fill_grey() +   
  theme(legend.position="none")
print(robbery)


# b. Assault ----

figdata <- Rmisc::summarySE(data, 
                            measurevar="assault", 
                            groupvars=c("lottery_winner"))

figdata$point <- sprintf(figdata$assault*100, fmt = '%#.2f')

# p-val ri
ITT_CACE_pri <- paste0(#"CACE: ", round(CACE_mechs[['Assault']]$coefficients[2],3)*100,"%p\n",
  "ITT: ", round(ri_mechs[['Assault']]$estimate,3)*100,"%p\n",
  "(p=", round(ri_mechs[['Assault']]$two_tailed_p_value,3),")")

assault <- ggplot(figdata, aes(x=factor(lottery_winner), 
                               y=100*assault, 
                               fill=factor(lottery_winner))) +
  geom_bar(stat="identity", width=0.45, alpha=0.9) +
  geom_text(aes(label = point), size = 4.0, vjust=-0.5, hjust = 0.4, nudge_x = 0) +
  xlab("") +
  ylab(" ") +
  ylim(-2, 27) +
  scale_x_discrete(breaks=c("0", "1"), labels=mylabels)+
  theme_minimal() +
  theme(axis.text.y=element_text(colour = "gray30", size=11), axis.ticks.y=element_blank()) +
  theme(axis.text.x=element_text(colour = "gray30", size=11), axis.ticks.x=element_line(colour="gray30")) + 
  geom_segment(aes(x=1, y=19, xend=1.2, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.2, y=21, xend=1.8, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.8, y=21, xend=2, yend=19), colour="darkgray") + 
  annotate("text", x=1.5, y=21, label=ITT_CACE_pri, size=4) +  
  scale_fill_grey() + 
  theme(legend.position="none")
print(assault)


# c. Discretionary ----

figdata <- Rmisc::summarySE(data, 
                            measurevar="discretionary", 
                            groupvars=c("lottery_winner"))

figdata$point <- sprintf(figdata$discretionary*100, fmt = '%#.2f')

# ri p-val
ITT_CACE_pri <- paste0(#"CACE: ", round(CACE_mechs[['Discretionary']]$coefficients[2],3)*100,"%p\n",
  "ITT: ", round(ri_mechs[['Discretionary']]$estimate,3)*100,"%p\n",
  "(p=", round(ri_mechs[['Discretionary']]$two_tailed_p_value,3),")")

discretionary <- ggplot(figdata, aes(x=factor(lottery_winner), 
                                     y=100*discretionary, 
                                     fill=factor(lottery_winner))) +
  geom_bar(stat="identity", width=0.45, alpha=0.9) +
  geom_text(aes(label = point), size = 4.0, vjust=-0.5, hjust = 0.4, nudge_x = 0) +
  xlab("") +
  ylab(" ") +
  ylim(-2, 27) +
  scale_x_discrete(breaks=c("0", "1"), labels=mylabels)+
  theme_minimal() +
  theme(axis.text.y=element_text(colour = "gray30", size=11), axis.ticks.y=element_blank()) +
  theme(axis.text.x=element_text(colour = "gray30", size=11), axis.ticks.x=element_line(colour="gray30")) + 
  geom_segment(aes(x=1, y=19, xend=1.2, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.2, y=21, xend=1.8, yend=21), colour="darkgray") +
  geom_segment(aes(x=1.8, y=21, xend=2, yend=19), colour="darkgray") + 
  annotate("text", x=1.5, y=21, label=ITT_CACE_pri, size=4) +   
  scale_fill_grey() + 
  theme(legend.position="none")
print(discretionary)

# d. Final figure ----
Figure_5<-ggarrange(robbery, assault, discretionary, 
                    ncol=3, nrow=1, labels = c("Robbery","Assault", "Discretionary"), 
                    font.label = list(size = 15))

pdf("Output/Figure5.pdf",height=5, width=10)
Figure_5
dev.off()
